Projet Chamois

1 Chargement des librairies


library(tidyverse)
library(corrplot)
library(lmerTest)
library(ade4)
library(splines)
library(car)
library(plotly)
library(DT)
library(Hmisc)
library(kableExtra)
library(knitr)

2 Import et description du jeu de données


2.1 Import des données

2.2 Description des données


2.2.1 Résumé des données

## 'data.frame':    1328 obs. of  7 variables:
##  $ id    : Factor w/ 217 levels "101","105","106",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ year  : int  1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 ...
##  $ fec   : int  1 1 1 1 1 1 1 0 0 0 ...
##  $ coh   : int  1995 1995 1995 1995 1995 1995 1995 1995 1995 1995 ...
##  $ anmark: int  1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 ...
##  $ pds   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ ydth  : int  2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
## cham 
## 
##  7  Variables      1328  Observations
## --------------------------------------------------------------------------------
## id 
##        n  missing distinct 
##     1328        0      217 
## 
## lowest : 101 105 106 107 108, highest: 82  87  9   93  R1 
## --------------------------------------------------------------------------------
## year 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1328        0       27    0.998     2006    6.831     1995     1997 
##      .25      .50      .75      .90      .95 
##     2001     2006     2010     2014     2015 
## 
## lowest : 1991 1992 1993 1994 1995, highest: 2013 2014 2015 2016 2017
## --------------------------------------------------------------------------------
## fec 
##        n  missing distinct     Info      Sum     Mean      Gmd 
##     1328        0        2    0.716      806   0.6069   0.4775 
## 
## --------------------------------------------------------------------------------
## coh 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1328        0       33    0.997     1996     7.75     1985     1987 
##      .25      .50      .75      .90      .95 
##     1991     1997     2001     2005     2007 
## 
## lowest : 1977 1978 1980 1982 1983, highest: 2007 2009 2010 2011 2014
## --------------------------------------------------------------------------------
## anmark 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1328        0       24    0.996     2002    6.288     1993     1994 
##      .25      .50      .75      .90      .95 
##     1998     2002     2006     2009     2011 
## 
## lowest : 1991 1992 1993 1994 1995, highest: 2010 2011 2012 2014 2015
## --------------------------------------------------------------------------------
## pds 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     1100      228       92    0.999    19.89     5.25     11.5     12.0 
##      .25      .50      .75      .90      .95 
##     16.9     21.1     23.3     25.0     26.0 
## 
## lowest :  7.8 10.5 11.0 11.1 11.3, highest: 26.5 26.8 27.0 28.3 28.4
## --------------------------------------------------------------------------------
## ydth 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      920      408       22    0.977     2006    4.908     1998     2000 
##      .25      .50      .75      .90      .95 
##     2003     2007     2008     2012     2014 
## 
## lowest : 1994 1996 1997 1998 1999, highest: 2012 2013 2014 2015 2016
## --------------------------------------------------------------------------------

2.2.2 Histogramme années d’observation

2.2.3 Histogramme individus suivis

2.2.4 Présentation des données

Le data frame est constitué de 7 variables et 1328 observations. Chaque observation correspond à l’information de fécondité associée à une femelle chamois et relative à une année donnée. Le jeu de données résume les suivis réalisés entre 1991 et 2017 sur 27 années. D’après l’histogramme des années d’observation, les années entre 2004 et 2006 sont les années pour lesquelles le nombre de chamois suivis a été le plus important (cf histogramme années d’observation). 217 femelles chamois ont été suivies au total. Le nombre d’années de suivi varie selon les femelles entre 1 et 16 années (cf histogramme individus suivis).

2.3 Prétraitement du jeu de données

2.3.1 Elimination des données aberrantes

Les chamois observés après leur mort ou avant leur naissance sont retirés du jeu de donnees.

cham <- cham %>% 
  filter(year<=ydth | is.na(cham$ydth)) %>%
  filter(year>=coh)

2.3.2 Création des variables âge (age) et longévité (long)

cham2 <- cham %>%
  summarise(cham, age= year-coh, long=ydth-coh, agemark=anmark-coh)

3 Question 1 : Lien fécondité annuelle et âge des femelles


3.1 Représentation graphique des données

3.1.1 Représentation par classe d’âge

3.1.2 Représentation sans grouper par classe d’âge

3.1.2.1 Utilisation de la fonction jitter

3.1.2.2 Utilisation de la fonction geom_count

Graphiquement, une augmentation de l’âge des chamois semble engendrer une diminution de la fécondité des chamois avec une inflexion de la courbe observée autour de 10 ans.

3.2 Analyse statistique du lien entre fécondité annuelle et l’âge des femelles

3.2.1 Modèles de régression lineaire generalisé avec effets aléatoires


3.2.1.1 Premier modèle testé

Le premier modèle appliqué est un modèle glm qui utilise la fonction de lien binomial afin de prendre en compte le fait que la variable réponse soit une variable binomiale. La variable “id” est désignée comme variable aléatoire pour prendre en compte le fait que les observations sont répetées sur les mêmes individus sur plusieurs années.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ age + (1 | id)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1761.9   1777.4   -877.9   1755.9     1314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6385 -1.1378  0.6862  0.7921  1.0431 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.2088   0.4569  
## Number of obs: 1317, groups:  id, 217
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.58495    0.15766   3.710 0.000207 ***
## age         -0.01789    0.01565  -1.143 0.253095    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## age -0.903
Chisq Df Pr(>Chisq)
age 1.30614 1 0.2530947

Interpretation des coefficients:

## [1] 1.336377
##                Estimate Std. Error   z value     Pr(>|z|)
## (Intercept)  0.58494620 0.15766366  3.710089 0.0002071865
## age         -0.01788695 0.01565097 -1.142865 0.2530947451
## [1] 0.9822721

L’AIC de ce modèle est de 1736. Avec ce modele, la dispersion calculée comme le ratio variance/df est de 1.3 donc il n’y a pas de surdispersion observee. On calcule en utilisant la fonction inverse de logit le coefficient qui permet d’exprimer la fécondité annuelle (en prenant en compte la fonction de lien) en fonction de l’âge. Il est 1.8047875% moins vraisemblable que les chamois aient un petit lorsque leur âge augmente d’un an mais la p value est de 0.25 donc l’impact de l’âge via ce premier modèle n’est pas significatif.

3.2.1.2 Second modèle testé

Un modèle quadratique est testé par la suite pour prendre en compte la tendance de la ligne de régression observée sur les graphiques qui présente une inflexion autour de 10 ans.

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00435823 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ age + I(age^2) + (1 | id)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1674.2   1694.9   -833.1   1666.2     1313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1557 -1.0113  0.5808  0.7269  2.8480 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.2407   0.4906  
## Number of obs: 1317, groups:  id, 217
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.076697   0.338088  -6.142 8.12e-10 ***
## age          0.642350   0.077437   8.295  < 2e-16 ***
## I(age^2)    -0.033970   0.003999  -8.495  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) age   
## age      -0.948       
## I(age^2)  0.875 -0.977
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00435823 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

La variable age est centrée normée car le modèle n’arrive pas à converger.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ age_scale + I(age_scale^2) + (1 | id)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1674.2   1694.9   -833.1   1666.2     1313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1557 -1.0113  0.5808  0.7269  2.8481 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.2407   0.4906  
## Number of obs: 1317, groups:  id, 217
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     0.95576    0.09371  10.199   <2e-16 ***
## age_scale       0.09431    0.06654   1.417    0.156    
## I(age_scale^2) -0.53218    0.06267  -8.492   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ag_scl
## age_scale    0.146       
## I(ag_scl^2) -0.654 -0.168

Interpretation des coefficients:

## [1] 1.26885
##                   Estimate Std. Error   z value     Pr(>|z|)
## (Intercept)     0.95576463 0.09371212 10.198944 2.004403e-24
## age_scale       0.09430895 0.06653992  1.417329 1.563868e-01
## I(age_scale^2) -0.53217960 0.06267034 -8.491730 2.035830e-17

L’AIC de ce modèle est de 1674. Avec ce modèle, la dispersion calculée est de 1.3 donc il n’y a pas de surdispersion observée. L’AIC de ce modèle 2 quadratique < l’AIC du modèle 1 linéaire donc le modèle quadratique est plus adapté comme attendu graphiquement. Une observation des coefficients associés aux termes âge et âge^2 indique que le terme “âge” n’est pas significatif dans la prédiction de la variable réponse alors que la p value associée au terme “âge^2” < 0.01. La fonction carré est donc testée.

3.2.1.3 Troisième modèle testé

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ I(age_scale^2) + (1 | id)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1674.2   1689.7   -834.1   1668.2     1314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1199 -1.0427  0.5836  0.7341  3.0118 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.2536   0.5036  
## Number of obs: 1317, groups:  id, 217
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     0.93801    0.09350   10.03   <2e-16 ***
## I(age_scale^2) -0.51897    0.06283   -8.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## I(ag_scl^2) -0.647

Interpretation des coefficients:

## [1] 1.269406
##                  Estimate Std. Error   z value     Pr(>|z|)
## (Intercept)     0.9380093 0.09350033 10.032150 1.100933e-23
## I(age_scale^2) -0.5189717 0.06282985 -8.259955 1.457278e-16
## [1] 0.5951322

L’AIC de ce modèle est de 1674. Avec ce modèle, la dispersion calculée est 1.3 donc il n’y a pas de surdispersion observée. Il est 1.8047875% moins vraisemblable que les chamois aient un petit lorsque leur âge au carré augmente d’une unité (p value < 0.05).

3.2.1.4 4ème modèle testé

On rajoute la variable “year” comme variable aléatoire pour prendre également en compte le fait que les mêmes individus ont été suivis sur les mêmes années.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ I(age_scale^2) + (1 | id) + (1 | year)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1640.9   1661.6   -816.4   1632.9     1313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2022 -0.8931  0.5182  0.6924  5.6258 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.2863   0.5351  
##  year   (Intercept) 0.3040   0.5514  
## Number of obs: 1317, groups:  id, 217; year, 27
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.01678    0.14911   6.819 9.17e-12 ***
## I(age_scale^2) -0.56448    0.06598  -8.555  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## I(ag_scl^2) -0.431

Interpretation des coefficients:

## [1] 1.243717

L’AIC de ce modèle est de 1641. Avec ce modèle, la dispersion calculée est de 1.2 donc il n’y a pas de surdispersion observée.

3.2.2 Conclusions

anova(glm1, glm1q, glm1c, glm1d)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
glm1 3 1761.863 1777.412 -877.9314 1755.863 NA NA NA
glm1c 3 1674.160 1689.710 -834.0802 1668.160 87.702339 0 NA
glm1q 4 1674.163 1694.895 -833.0813 1666.163 1.997873 1 0.1575201
glm1d 4 1640.892 1661.624 -816.4458 1632.892 33.270918 0 NA

Le modèle glm1d présente le plus faible AIC. La variable “âge” présente un effet significatif sur la fécondité annuelle via ce modèle ce qui n’est pas surprenant car graphiquement la ligne de régression présentait une courbe avec une diminution de la fécondité pour des âges élevés.

4 Variation de la fécondité annuelle en fonction du temps


4.1 Représentation graphique des données

4.1.1 Représentation graphique par année

4.1.2 Représentation graphique sans grouper par annee

Graphiquement, la variable âge ne semble pas avoir d’effet sur la fécondité annuelle des chamois malgré le fait que l’âge moyen de la population augmente sensiblement avec les années.

4.2 Analyse statistique du lien entre fécondité annuelle et années

4.2.1 Modèles de régression lineaire generalisé avec effets aléatoires


4.2.1.1 Premier modèle

Le premier modèle appliqué est un modèle glm qui utilise la fonction de lien binomial afin de prendre en compte le fait que la variable réponse soit une variable binomiale. La variable “id” est désignée comme variable aléatoire pour prendre en compte le fait que les observations sont répetées sur les mêmes individus sur plusieurs années.

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.28972 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ year_scale + (1 | id)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1763.1   1778.7   -878.6   1757.1     1314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6594 -1.1429  0.6860  0.7915  1.0701 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.2138   0.4623  
## Number of obs: 1317, groups:  id, 217
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.42246    0.06787   6.225 4.83e-10 ***
## year_scale  -0.01480    0.06540  -0.226    0.821    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## year_scale 0.002

Interpretation des coefficients:

## [1] 1.337139

L’AIC de ce modèle est de 1763. Avec ce modèle, la dispersion calculée comme le ratio variance/df est de 1.3 donc il n’y a pas de surdispersion observée. D’après la p-value > 0.82, il n’y a pas d’effets significatifs de la variable “year” sur la fécondité annuelle comme supposé préalablement par les représentations graphiques. La variable “âge” n’ayant pas d’effets significatifs sur la fécondité annuelle des chamois, il n’y a pas d’utilité à appliquer des modèles avec effets additifs ou multiplicatifs en utilisant les variables âge et année.

5 Lien entre fecondite totale et longevite des animaux


5.1 Représentation graphique des données

5.1.1 Représentation sans prendre en compte le nombre d’annees de suivi

cham_id = cham2 %>% 
  group_by(id) %>% 
  dplyr::summarise(feconditetotale= sum(fec), long=long, pds=pds, coh=coh, anneetot=(ydth-min(year)+1), MinAn=min(year), MaxAn=max(year), AgeMark=min(age)) %>%
  unique()%>%
  drop_na(long)
plot9 <-ggplot(cham_id, aes(x=long, y=feconditetotale)) +
    geom_count() +
    labs(title = "Fecondite totale en fonction de la longevite",x="Longevite", y="Fecondite totale") +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_smooth(); plot9

plot10 <- ggplot(cham_id, aes(x=long, y=feconditetotale)) +
    geom_jitter(width = 0.25, height = 0.25)+
     labs(title = "Fecondite totale en fonction de la longevite",x="Longevite", y="Fecondite totale") +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_smooth(); plot10

### Prise en compte du biais apporté par le nombre d’années de suivi

Tous les chamois n’ont pas été marqués au même age et donc n’ont pas été suivis le même nombre d’annee comme le montre l’histogramme ci-dessous.

plot11 <- ggplot(cham_id, aes(x=anneetot)) +
    geom_histogram(bins = 30)+
     labs(title = "Nombre annees de suivi des chamois",x="Nombre d'annees") +
    theme(plot.title = element_text(hjust = 0.5)); plot11

plot(jitter(cham_id$long), jitter(cham_id$anneetot))

Le nombre d’annees de suivi n’est donc pas égal à la longévité des individus comme illustré par le graphique ci-dessous.

plot5 <- ggplot(data = cham_id, aes(x = anneetot,y=long))+
     labs(title = "Correlation entre le nombre d'annees de suivi et la longevite",x="Annees de suivi", y="Longevite") +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_point()+
    geom_smooth(); plot5

On vérifie dans un premier temps, l’impact du nombre d’annees de suivi sur la fecondite totale des chamois

plot10bis <- ggplot(cham_id, aes(x=anneetot, y=feconditetotale)) +
    geom_point(stat="identity") +
    labs(title = "Fecondite totale en fonction du nombre d'annee de suivi",x="Nombre annees de suivi", y="Fecondite totale") +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_smooth(); plot10bis

Comme attendu, la fecondite totale augmente avec le nombre d’annees de suivi. Or il est difficile de savoir si cette augmentation est du au fait que l’individu a été marqué précocement ou si l’individu a vécu plus longtemps.

Pour pouvoir répondre à la question initiale, qui consiste à vérifier s’il y a un lien entre la fécondité annuelle et la longévité, il faut pouvoir comparer des individus suivis sur le même nombre d’année et si possible sur le maximum d’années possibles.

plot11 <- ggplot(cham_id, aes(x=AgeMark)) +
    geom_histogram(bins=50)+
     labs(title = "Age de marquage des individus",x="Age") +
    theme(plot.title = element_text(hjust = 0.5)); plot11

En s’appuyant sur le graphique ci-dessous, l’analyse statistique va être réalisé sur un échantillon de la population suivi à savoir tous les individus marqués à 3 ans qui constitue un échantillon > 30 individus et peut permettre d’étudier l’impact de la longévité sur la fécondité totale.

On sait que la variable “année” n’a pas d’impact sur la fécondité annuelle des chamois donc le faut que les chamois aient été suivis pendant des périodes différentes n’a pas d’impact. Il faut par contre vérifier au préalable que le fait de sélectionner la fécondité de la population des individus marqués à 3 ans est représentative de la population totale de chamois.

r <- tapply(cham2$fec, cham2$agemark, mean)
r
##         0         1         2         3         4         5         6         7 
## 0.5784314 0.6170213 0.6818182 0.6699029 0.6346154 0.5725806 0.6608696 0.6477273 
##         8         9        10        11        12        13        14        15 
## 0.6071429 0.5740741 0.5573770 0.6065574 0.5483871 0.5454545 0.5000000 1.0000000 
##        16        17        21 
## 0.2000000 0.0000000 0.0000000
Age3 <- cham2[cham2$agemark==3,]
Ageautre <- cham2[cham2$agemark!=3,]
mean(Age3$fec)
## [1] 0.6699029
mean(Ageautre$fec)
## [1] 0.6004942
#Voir avec Karim quel test

5.1.2 Création du sous-échantillon pour répondre à la question

cham3 <- cham_id %>% 
  filter(AgeMark==3)

5.1.3 Représentation graphique

plot9 <-ggplot(cham3, aes(x=long, y=feconditetotale)) +
    geom_point() +
    labs(title = "Fecondite totale en fonction de la longevite",x="Longevite", y="Fecondite totale") +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_smooth(); plot9

5.2 Analyse statistique du lien entre la fécondité annuelle et la longevite

5.2.1 Tests de modèles de régression lineaire generalise avec effets aléatoires


5.2.1.1 First

Plusieurs modèles ont été testés => voir avec Karim

lm1 <-lm(feconditetotale~long, data = cham3)
summary(lm1)
## 
## Call:
## lm(formula = feconditetotale ~ long, data = cham3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9408 -1.1635  0.2819  1.0436  3.5826 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.52025    0.64088  -3.932 0.000275 ***
## long         0.74610    0.06406  11.648 1.86e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.759 on 47 degrees of freedom
## Multiple R-squared:  0.7427, Adjusted R-squared:  0.7372 
## F-statistic: 135.7 on 1 and 47 DF,  p-value: 1.863e-15
plot(lm1)

glmp <- glm(data=cham3, feconditetotale~long, family="poisson")
summary(glmp)
## 
## Call:
## glm(formula = feconditetotale ~ long, family = "poisson", data = cham3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3726  -0.5124  -0.1843   0.6590   1.3798  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.33487    0.22434  -1.493    0.136    
## long         0.17122    0.01803   9.494   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 133.056  on 48  degrees of freedom
## Residual deviance:  37.754  on 47  degrees of freedom
## AIC: 186.32
## 
## Number of Fisher Scoring iterations: 4
lm2 <- lm(feconditetotale~log(long), data = cham3)
summary(lm2)
## 
## Call:
## lm(formula = feconditetotale ~ log(long), data = cham3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3986 -1.5727  0.0741  1.2620  5.0012 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.3397     1.3263  -5.534 1.36e-06 ***
## log(long)     5.5322     0.6116   9.045 7.38e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.095 on 47 degrees of freedom
## Multiple R-squared:  0.6351, Adjusted R-squared:  0.6274 
## F-statistic: 81.81 on 1 and 47 DF,  p-value: 7.377e-12
plot(lm2)

6 Lien entre fecondite annuelle et longevite des animaux


6.1 Représentation graphique des données

plot12 <- ggplot(cham2, aes(x=long, y=fec)) +
    geom_count() + 
    labs(title = "Fécondité annuelle en fonction de la longevite",x="Longevite", y="Fécondité annuelle")+
    geom_smooth()+
    theme(plot.title = element_text(hjust = 0.5)); plot12

plot13 <- ggplot(cham2, aes(x=long, y=fec)) +
    geom_jitter(width = 0.55, height=0) + 
    labs(title = "Fécondité annuelle en fonction de la longevite",x="Longevite", y="Fécondité annuelle")+
    geom_smooth()+
    theme(plot.title = element_text(hjust = 0.5)); plot13

L’allure concave des lignes de régression illustre une augmentation de la fécondité annuelle avec la longévité jusqu’à atteindre un maximum autour de 13-14 ans puis une diminution de la fécondité annuelle.

6.2 Analyse statistique du lien entre fecondite annuelle et longevite des femelles

6.2.1 Modèles de régression lineaire generalise avec effets aléatoires


6.2.1.1 First

Le premier modele applique est un modele glm qui utilise la fonction de lien binomiale afin de prendre en compte le fait que la variable réponse soit une variable binomiale. La variable “id” est désignée comme variable aléatoire pour prendre en compte le fait que les observations sont repetees sur les mêmes individus sur plusieurs annees.

glm11 <- glmer(fec ~ long + (1| id),data=cham2, family = binomial)
long_scale <- scale(x = cham2$long,center = TRUE,scale = TRUE)
glm11 <- glmer(fec ~ long_scale + (1| id),data=cham2, family = binomial)
summary(glm11)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ long_scale + (1 | id)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1218.0   1232.5   -606.0   1212.0      906 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7916 -1.0634  0.6517  0.7850  1.2047 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.3199   0.5656  
## Number of obs: 909, groups:  id, 163
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.38681    0.08681   4.456 8.36e-06 ***
## long_scale   0.04500    0.08512   0.529    0.597    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## long_scale 0.095
Anova(glm11)
Chisq Df Pr(>Chisq)
long_scale 0.2793972 1 0.5970966

Interpretation des coefficients:

dispersion_glm1 <- 1212/906;dispersion_glm1 #Surdispersion = deviance/df
## [1] 1.337748

L’AIC de ce modèle = 1218. Avec ce modele, la dispersion calculée comme le ratio variance/df = 1.3 donc il n’y a pas de surdispersion observee. Avec ce modele, la p value associé à l’impact de la varibale “longevite” sur la fecondite annuelle est de 0.6 donc l’effet de la longevite sur la variable reponse n’est pas significatif.

6.2.1.2 Second

On applique un modèle quadratique pour prendre en compte la tendance la ligne de régression observee sur les graphiques qui présente une inflexion autour de 13-14 ans.

glm11 <- glmer(fec ~ long_scale + I(long_scale^2)+ (1| id),data=cham2, family = binomial)
summary(glm11)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ long_scale + I(long_scale^2) + (1 | id)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1209.3   1228.6   -600.7   1201.3      905 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7429 -1.0969  0.6542  0.7821  1.4323 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.2492   0.4992  
## Number of obs: 909, groups:  id, 163
## 
## Fixed effects:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.60021    0.10574   5.676 1.38e-08 ***
## long_scale       0.03059    0.08331   0.367  0.71352    
## I(long_scale^2) -0.20676    0.06313  -3.275  0.00106 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lng_sc
## long_scale   0.026       
## I(lng_sc^2) -0.612  0.065
Anova(glm11)
Chisq Df Pr(>Chisq)
long_scale 0.1347913 1 0.7135150
I(long_scale^2) 10.7277777 1 0.0010554

Interpretation des coefficients:

dispersion_glm1 <- 1201/905;dispersion_glm1 #Surdispersion = deviance/df
## [1] 1.327072

L’AIC de ce modèle glm1q = 1209. Avec le modele glm1q, la dispersion calculée = 1.3 donc il n’y a pas de surdispersion observee. L’AIC de ce modèle glm1q < l’AIC du modèle glm1 donc le modèle glm1q est plus adapté. Une observation des coefficients indique que le terme “x” n’est pas significatif dans la prediction de la variable reponse alors que la p value associée au terme “x2” < 0.01. La fonction carré est donc testée.

6.2.1.3 Third

glm11 <- glmer(fec ~ I(long_scale^2)+ (1| id),data=cham2, family = binomial)
summary(glm11)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ I(long_scale^2) + (1 | id)
##    Data: cham2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1207.5   1221.9   -600.7   1201.5      906 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7388 -1.0973  0.6529  0.7774  1.3841 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.2565   0.5064  
## Number of obs: 909, groups:  id, 163
## 
## Fixed effects:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.59931    0.10618   5.644 1.66e-08 ***
## I(long_scale^2) -0.20842    0.06328  -3.294 0.000989 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## I(lng_sc^2) -0.616
Anova(glm11)
Chisq Df Pr(>Chisq)
I(long_scale^2) 10.84755 1 0.0009893

Interpretation des coefficients:

dispersion_glm1c <- 1201/906;dispersion_glm1c #Surdispersion = deviance/df
## [1] 1.325607
summary(glm11)$coefficients
##                   Estimate Std. Error   z value     Pr(>|z|)
## (Intercept)      0.5993126 0.10618061  5.644275 1.658785e-08
## I(long_scale^2) -0.2084191 0.06328076 -3.293563 9.892629e-04
exp(summary(glm11)$coefficients[[2]])
## [1] 0.8118667
Coeffb <- ((1/exp(summary(glm11)$coefficients[2]))-1)*100

L’AIC de ce modèle glm1c = 1674. Avec le modele glm1c, la dispersion calculée = 1.3 donc il n’y a pas de surdispersion observee. Il est 23.1729325% moins vraisemblable que les chamois aient un petit lorsque leur longevite au carré augmente d’une unité (p value < 0.05).

Segmentation => voir avec Karim

cham_long1 <- cham2[cham2$long<15,]
cham_long2 <- cham2[cham2$long>=15,]

plot16 <- ggplot(cham_long1, aes(x=long, y=fec)) +
    geom_count() + 
    labs(title = "Fécondité en fonction de la longevite",x="Longevite", y="Fécondité")+
    geom_smooth()+
    theme(plot.title = element_text(hjust = 0.5)); plot16
## Warning: Removed 408 rows containing non-finite values (`stat_sum()`).
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 408 rows containing non-finite values (`stat_smooth()`).

glm12 <- glmer(fec ~ long + (1| id),data=cham_long1, family = binomial)
summary(glm12)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ long + (1 | id)
##    Data: cham_long1
## 
##      AIC      BIC   logLik deviance df.resid 
##    789.0    802.2   -391.5    783.0      585 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6864 -1.0755  0.6737  0.8230  1.1503 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.1568   0.3959  
## Number of obs: 588, groups:  id, 120
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.66433    0.39918  -1.664  0.09606 . 
## long         0.09835    0.03665   2.684  0.00728 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## long -0.971
Anova(glm12)
Chisq Df Pr(>Chisq)
long 7.201696 1 0.0072835
plot17 <- ggplot(cham_long2, aes(x=long, y=fec)) +
    geom_count() + 
    labs(title = "Fécondité en fonction de la longevite",x="Longevite", y="Fécondité")+
    geom_smooth()+
    theme(plot.title = element_text(hjust = 0.5)); plot17
## Warning: Removed 408 rows containing non-finite values (`stat_sum()`).
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 408 rows containing non-finite values (`stat_smooth()`).

glm13 <- glmer(fec ~ long + (1| id),data=cham_long2, family = binomial)
summary(glm13)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: fec ~ long + (1 | id)
##    Data: cham_long2
## 
##      AIC      BIC   logLik deviance df.resid 
##    418.9    430.3   -206.5    412.9      318 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9648 -1.0836  0.5871  0.7880  1.2014 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.3949   0.6284  
## Number of obs: 321, groups:  id, 43
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  4.35885    1.43312   3.041  0.00235 **
## long        -0.23441    0.08423  -2.783  0.00539 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## long -0.994
Anova(glm13)
Chisq Df Pr(>Chisq)
long 7.744764 1 0.0053868

7 Lien entre fecondite totale et poids


7.1 Représentation graphique des données

7.1.1 Vérification de la comparabilite des poids selon les ages de capture

cham_pds <- cham_id %>% 
  drop_na(pds)
plot18 <- ggplot(cham_pds, aes(x=AgeMark, y=pds)) + 
    geom_point(color="blue", fill=rgb(0.1,0.4,0.5,0.7)) + 
    labs(title = "Poids selon l'age de capture",x="Age de capture", y="Poids mesuré")+
    theme(plot.title = element_text(hjust = 0.5))+
    geom_smooth();plot18

#Fiare une boucle for pour exclure outlier
(boxplot(cham_pds$pds))

## $stats
##      [,1]
## [1,] 16.0
## [2,] 20.0
## [3,] 22.0
## [4,] 23.5
## [5,] 27.0
## 
## $n
## [1] 131
## 
## $conf
##          [,1]
## [1,] 21.51684
## [2,] 22.48316
## 
## $out
##  [1] 14.3 12.8  7.8 14.5 12.5 14.0 14.2 11.9 12.5 14.0 11.1 11.3 11.1 14.5 11.5
## [16] 10.5
## 
## $group
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $names
## [1] "1"
out <- boxplot.stats(cham_pds$pds)$out
max(out)
## [1] 14.5
cham_pds <- cham_pds%>%
  filter(pds>max(out))

hist(cham_pds$AgeMark)

plot(cham_pds$pds, cham_pds$feconditetotale)

plot(cham_pds$AgeMark, cham_pds$feconditetotale)

Reflechir à la section des donnes avant analyse stat pour fecondite totale

7.2 Analyse statistique du lien entre fecondite annuelle/longevite et poids des femelles

7.2.1 Modèles de régression lineaire generalise avec effets aléatoires


7.2.1.1 First

lm14 <- lm(feconditetotale ~ pds,data=cham_pds)
summary(lm14)
## 
## Call:
## lm(formula = feconditetotale ~ pds, data = cham_pds)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.205 -2.086 -0.979  1.790  9.910 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.28615    2.26237   1.011    0.314
## pds          0.03828    0.10097   0.379    0.705
## 
## Residual standard error: 2.526 on 113 degrees of freedom
## Multiple R-squared:  0.00127,    Adjusted R-squared:  -0.007568 
## F-statistic: 0.1437 on 1 and 113 DF,  p-value: 0.7053
lm15 <- lm(long ~ pds,data=cham_pds)
summary(lm15)
## 
## Call:
## lm(formula = long ~ pds, data = cham_pds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.1325 -2.3945  0.0646  2.3182  7.8675 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.4877     3.0358  -0.490    0.625    
## pds           0.6010     0.1355   4.435 2.14e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.389 on 113 degrees of freedom
## Multiple R-squared:  0.1483, Adjusted R-squared:  0.1407 
## F-statistic: 19.67 on 1 and 113 DF,  p-value: 2.143e-05

7.2.2 Représentation graphique de la fecondite totale en fonction du poids

plot19 <- ggplot(cham_pds, aes(x=pds, y=feconditetotale)) + 
    geom_point(color="blue", fill=rgb(0.1,0.4,0.5,0.7)) + 
    geom_smooth();plot19

7.2.3 Représentation graphique de la longevite en fonction du poids

plot20 <- ggplot(cham_pds, aes(x=pds, y=long)) + 
    geom_point(color="blue", fill=rgb(0.1,0.4,0.5,0.7)) + 
    geom_smooth();plot20

Longevite significatif sur le poids mais pas de la fecondite totale.